Rotation induced superfluid-normal phase separation in trapped Fermi gases 
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We use the Bogoliubov-de Gennes formalism to analyze the effects of rotation on the ground state phases of 
harmonically trapped Fermi gases, under the assumption that quantized vortices are not excited. We find that 
the rotation breaks Cooper pairs that are located near the trap edge, and that this leads to a phase separation 
between the nonrotating superfluid (fully paired) atoms located around the trap center and the rigidly rotating 
normal (nonpaired) atoms located towards the trap edge, with a coexistence (partially paired) region in between. 
Furthermore, we show that the superfluid phase that occurs in the coexistence region is characterized by a 
gapless excitation spectrum, and that it is distinct from the gapped phase that occurs near the trap center. 
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With the ultimate success of techniques for trapping and 
cooling atomic gases developed and improved gradually since 
the 1980s, atomic Fermi gases have emerged as unique testing 
grounds for many theories of exotic matter in nature, allowing 
for the creation of complex yet very accessible and control- 
lable many-body quantum systems. For instance, evolution 
from the weakly attracting Bardeen-Cooper-Schrieffer (BCS) 
limit to the weakly repulsive molecular Bose-Einstein conden- 
sation (BEC) have been observed in a series of remarkable 
experiments. The ground state phases of such atomic Fermi 
gases have since been the subject of intense theoretical and 
experimental research worldwide iH, . 

To verify the superfluid ground state of atomic Fermi gases, 
it is essential to analyze their response to rotation. For in- 
stance, it has been theoretically predicted that sufficiently fast 
rotation of superfluid Fermi gases excites quantized vortices 
in the form of hexagonal vortex lattices |3, 4, 5]. Such vor- 
tex lattices have recently been observed for both population 
balanced |6] and imbalanced |7] systems. These vortex ex- 
periments have not only complemented previously found sig- 
natures, but also provided the ultimate evidence to support the 
superfluid nature of the ground state. 

Recently, the effects of rotation on the ground state phases 
of harmonically trapped Fermi gases has been theoretically 
studied at unitarity, under the assumption that quantized vor- 
tices are not excited (see also Ref. (HI]). It has been argued 
that the rotation causes a complete phase separation between a 
nonrotating superfluid core and a rigidly rotating normal gas, 
with a discontinuous density at the interface. Subsequently, 
this problem has been analyzed using the mean-field BCS 
framework within the semi-classical local density approxima- 
tion (LDA) ilOll . In addition to the superfluid and normal re- 
gions, a coexistence region is found, the possibility of which 
is not considered in Ref. [8]. 

In this manuscrript, we go beyond the semi-classical LDA 
and develop a quantum mechanical Bogoliubov-de Gennes 
(BdG) formalism to analyze the effects of rotation on the 
ground state phases of harmonically trapped Fermi gases. We 
discuss both population balanced and imbalanced mixtures 
throughout the BCS-BEC evolution. Our main results are in 
qualitative agreement with those of Ref. [10]. We find that 
the rotation breaks Cooper pairs that are located near the trap 



edge, and that this leads to a phase separation between the 
nonrotating superfluid (fully paired) atoms located around the 
trap center and the rigidly rotating normal (nonpaired) atoms 
located towards the trap edge, with a coexisting (partially 
paired) region in between. This leads to a continuous density 
and superfluid order parameter as a function of radial distance. 
Furthermore, we show that the superfluid phase that occurs in 
the coexistence region is characterized by a gapless excitation 
spectrum, and that it is distinct from the gapped phase that 
occurs near the trap center 

We obtain these results by solving the mean-field BdG 
equations in the rotating frame (in units of ^ = fcs = 1) 

IC^{r)-nC, A(r) 1 [ ..^(r) 1 _ [ ..^(r) 

A*(r) -ICl{r) + nCt\ [v^{r) \ " [ i;,(r) 

(1) 

where ICa{r) = -V^/(2M) - fJ-air), ^ is the rotation fre- 
quency around the z-axis of the trapping potential, Cz is 
the z-component of the angular momentum operator, and the 
off-diagonal self-consistency field A(r) — g(?/'|(r)'(/;j(r)) is 
the local superfluid order parameter. Here, cr = {f, j} la- 
bels the trapped hyperfine states, M is the mass, /i(j(r) = 
fic — V{r) is the local chemical potential, fi^ is the global 
chemical potential, V{r) = Muj'^r'^ /2 is the trapping po- 
tential which is assumed to be spherically symmetric, u) is 
the trapping frequency, g > is the strength of the zero- 
ranged attractive interactions between t and | atoms, and 
(...) is a thermal average. The quasiparticle wavefunctions 
u,-i (r) and (r) are related to the particle annihilation opera- 
tor ■(/'o-(r) via the Bogoliubov-Valatin transformation ■0cr (r) = 
E^K,a(r)7,,,a - s^-y; ,,(r)7^ _^], where 7+ ,^ and 7,,^^ are 
the quasiparticle creation and annihilation operators, respec- 
tively, and S| ~ +1 and si — —1. Since the BdG equations 



are invariant under the transformation 



(r), 



""'f.i (i") ^ I i"^) ™d ~^ ~*^').T' i'- sufficient to solve 
only for u^{v) = (r), w,,(r) = u,,,x(r) and e,, = e^j as 
long as we keep all of the solutions with positive and negative 
eigenvalues. 

We assume A(r) = ^gJ2ri''^vi^)'''rii^)fi^v) ''^^^ with- 
out loosing generality, where f{x) ~ l/[exp(a;/T) + 1] 
is the Fermi function and T is the temperature. Further- 
more, we can relate g to the two-body scattering length ap 
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via l/g = -M/{ATTaF) + Mkc{r)/ {2tt'^) where fc2(r) = 
2M[€c + n{r)] and fi{r) = [/i|(r) + ^^(r)]/2. Here, is 
the energy cutoff to be specified below, and our results de- 
pend weakly on the particular value of ec provided that it 
is chosen sufficiently high. The order parameter equation 
has to be solved self-consistently with the number equation 
No- = J drna-{r), where no-(r) = {ijjl.{r)ijja-{r)) is the local 
density of fermions, leading to n^{r) = 

and nx(r) = EJ?^r,(r)P/(-e»;)- 

Next, we expand Mr,(r) and v,i{r) in the complete 
basis of the harmonic trapping potential eigenfunctions 
given by /C^(r)(/)„,f,„(r) = i^^g(j)n,e,m{r), where = 
uj(2n + £ + 3/2) — /io- is the eigenvalue and 0„,^,m(r) = 
Rn,e{i^)ye,m{0T, fr) is the eigenfunction. Here, n is 
the radial quantum number, and £ and m are the or- 
bital angular momentum and its projection, respectively. 
The angular part Yi,rn{Sr, fr) is a spherical harmonic and 
the radial part is Rn,i{r) = y^{Mujf/'^[n\/{n + £ + 
l/2)!]i/2e-*^'/2f'^Ll+^/^(f2), where f = y/MlJr is dimen- 
sionless and (x) is an associated Laguerre polynomial. 

We choose Cz = —id/dipr and 77 = {£, to,7}, lead- 
ing to Uf,™,T,(r) = Z^n C£,m,7,„(/)„,f,m(r) andw^^„,^(r) = 
X^n di^m,^ ,n<t>n,t,m{'^) ■ This expansion rcduces the BdG equa- 
tions to a 2{ni + 1) X 2{ni + 1) matrix eigenvalue problem 
for a given {£, m} state 
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Here, rii — {ric — £)/2 is the maximal radial quantum num- 
ber and Tic is the radial quantum number cutoff, such that 
we include only the single particle states with aj(2n + £ + 
3/2) < ec — LLj{nc + 3/2). In Eq. |2] the diagonal ma- 
trix element is K^'" = £(5n,n' where 6ij is the Kro- 

necker delta, and the off-diagonal matrix element is A"'" « 
/ r'^ dr A{r) Rn.i{r) Rn' ,e{r) . Although A(r) becomes axially 
symmetric when 51 7^ 0, it is convenient to define A(r) = 
/ dl7rA(r)/(47r) leading to 



A(r) 



_3_ 

47r 



E 

£,m,7,n,n' 



^lm,7,„('^)^im,7,«'(^)/(e^,™,7)>(3) 



where we introduced rI ^ (r) 



m,7 



,nRnA^) and 
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J (ir2r"'<T(r)/(47r) leading to 



Rn.e{r). Similarly, we define ncr{r) 



1 



■^,m,7,n,n' 



Lastly, Nfj reduces to = J2 



i,m,^,n ^,m,7,n 



/(Q,m.7) and 



Nl - E.,™,7,„4™,7,»/(-^^^".7)- When n ^ 0, the 
eigenfunction coefficients ce^m,-y,n and di^m.-y,n and the eigen- 
values 7 become independent of m, and Eqs. (|2]l, (|3]i 



and (|4ll reduce to the usual ones, see e.g. Ref. lilllillilH. 

Therefore, due to the coupling between different m states, the 
rotating case is numerically much more involved compared to 
the nonrotating case. 

In addition to ricr(r), we want to calculate the local den- 
sity of normal (rotating) fermions na-,N{r). For this pur- 
pose, we use the local current density Jcr(r) — {Jc{v)), 
where J^{y) ^ [l/(2A/i)][V^J(r)VV'a(r) - H.c] is the 
quantum mechanical probability current operator and H.c. 
is the Hermitian conjugate. This leads to J|(r) — 
[l/(2Mz)]E,KWV«r,(r)/(e,) - H.c] and 3^{v) = 
[l/(2Afi)]E^K(r)Vz;;(r)/(-e^) - H.c]. The current, 
similar to the classical case, can be written as Jcr(r) = 
?icr_jv(r)v(r), where 7icr,Af(i") is the local density and v(r) = 
riz X r is the local velocity of normal fermions correspond- 
ing to a rigid-body rotation. Since the normal fermions are 
expelled towards the trap edge, we approximate ncr,Ar(r) — 
J dflrUcr^N (r) / (47r) as 



na,N{r) 



^,m,7,n,n 



"^-Rf,™,7,„(^■) 



^?',m,7,n'W/('SaQ,m,7), (5) 



such that Jcr(r) = / df^r Jcr(r)/(47r) ~ f7rricr,Ar(r). 

Having discussed the BdG formalism, next, we analyze the 
ground state (T = 0) phases for both population balanced 
(TV-f — or /i| = /i|) and imbalanced (iV| ^ or 
/i| 7^ ^i) Fermi gases. This is achieved by solving the BdG 
equations (|2]i, ^ and (|4]l self-consistently as a function of 
the dimensionless parameter l/{kpaF) where /cf is the Fermi 
momentum defined via the Fermi energy 



ep = Lo{np + 3/2) = ^ - \mu:\1 « ^(37V) 



1/3 



(6) 



and N ^ N^^ + ^ {np + \){np + 2){nF + 3)/3. Here, 
np and r^? are the corresponding Fermi level and Thomas- 
Fermi radius, respectively. In our numerical calculations, we 
choose np = lb and ric = 65, which corresponds to a total 
of = 1632 fermions and €c ~ ^ep, respectively. Here, it 
is important to emphasize that we do not expect any qualita- 
tive change in our results with higher values of np and/or nc, 
except for minor quantitative variations. 

In Fig. [T] we consider a weakly interacting Fermi gas on 
the BCS side with l/{kFap) = —0.5, and show n„{r) and 
A(r) for nonrotating {Vt = 0) and rotating (fi = 0.5w) cases. 
When Ni = N^, we find that A(r) depletes everywhere in- 
side the trap and especially around the trap edge. This is be- 
cause it is easier to break the Cooper pairs that are located 
towards the trap edge in comparison to the ones that occupy 
the center (see below). Our result is quantitatively different 
from the LDA one where A(r) depletes only around the trap 
edge 1 10]. When — SN^, in addition to such an effect, 
the spatial modulation of A(r) disappears in the rotating case 
as shown in Fig.fljb). The depletion of A(7-) leads to a de- 
crease (increase) in n^ {r) around the trap center (edge) due 
to the centrifugal force caused by the rotation. It also leads 
to a phase separation between the nonrotating fully paired su- 
perfluid (FPS) atoms located around the trap center and the 
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(a): 1/(kp:aF) =-0.5 and Nt = N|, 
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FIG. 1: We compare the density ncr{r) and the superfluid order pa- 
rameter A(r) for nonrotating (SI = 0) and rotating (Q, — 0.5lj) 
Fermi gases. Here, l/(fcFfflF) = —0.5, and 7V| = 7V| in (a) and 
Nt; = 3A^j in (b). We also show the density nc.N{r) of normal 
fermions (squared-dotted line) for the rotating case. 



(a):1/(kpap) = 0.0 and = N| 




FIG. 2: We compare the density no-(r) and the superfluid order pa- 
rameter A(r) for nonrotating (SI = 0) and rotating (S7 — 0.7uj) 
Fermi gases. Here, l/{kFaF) = 0.0, and Ni = Ni in (a) and 
A^l = 3Ni in (b). We also show the density na,N{r) of normal 
fermions (squared-dotted line) for the rotating case. 



rigidly rotating normal (nonpaired) ones located towards the 
trap edge, with a coexistence region in between, which is in 
agreement with the LDA result (lu\. We characterize the FPS, 
coexistence and normal phases by na-{r) ^ n^^jsilr) — 0, 
ntriT) > ncr^iq{r) 7^ and ncr{r) — ncr^iq{r) ^ 0, respec- 
tively. Notice that A(r) decreases smoothly as a function of 
r from the FPS to the normal region where it vanishes, which 
is in contrast with the LDA result where A(r) is nonanalytic 
at the transition [10]. Since A(r) is finite in the coexistence 
region, this region corresponds to a partially paired superfluid 
(FPS), and it occupies a larger region compared to the LDA 
resuls [|10.1 . In addition, the trap center becomes a FPS for 
17 > 0.5u; when l/(fcFaF) = -0.5. 

In Fig. 12] we consider a strongly interacting Fermi gas at 
unitarity with l/(fci?ai?) = 0, and show na{r) and A(r) for 
nonrotating (17 = 0) and rotating (f7 = O.Toj) cases. The main 
effects of rotation are qualitatively similar to the weakly in- 
teracting case. However, since the Cooper pairs become more 
strongly bound as a function of the interaction strength, it re- 
quires much faster Q, to break them. For instance, at unitar- 
ity, the entire superfluid is robust for Q, < Q.3u!, and the trap 
center stays as an FPS even for ft ^ uj (not shown). There- 
fore, the effects of rotation become weaker as the interaction 
strength increases, and both the PPS and the normal regions 
eventually disappear in the molecular limit (not shown), i.e. 
rotation can not break any Cooper pair in the molecular limit. 

Another important observable is the local angular momen- 
tum defined by L^.a-i^^) — {i''l{''^)^z'^ry{'^)) ■ This leads to 
Lz^ir) = -iY,^u*^{-r)du^{v)/difrf{f-n) and Lz,i{v) = 
— i;^(r)9w*(r)/9(pr/(— e^). Since the superfluid atoms 



do not carry angular momentum, L^.cri'^) is directly re- 
lated to the local density of normal fermions via Lz.a-{r) = 
J dflrLz,a{r)/i'in) w Mflr'^n„,N{r)- Therefore, Lz,„{r) 
can be easily extracted from Figs. [1] and |2] For complete- 
ness, the total angular momentum ^ — J drLz^ai^) be- 
comes L^,t = E^,m,7,„,„' "^4"^,7,r^'/(ef,™,7) and L^,x = 

- T,Lm,j.n,n' "^^l^,-^ J {-u^m,^) , and It increases with in- 
creasing Vl. In the weakly attracting limit, L^.a becomes 
its rigid-body value when increases high enough so that 
A(r) — > everywhere (not shown). 

The microscopic mechanism responsible for the pair break- 
ing effects can be understood analytically within the semi- 
classical LDA, i.e. each component of the Fermi gas is con- 
sidered as locally homogenous at each position r with a lo- 
cal chemical potential /icr(r). In this approximation, the local 
quasiparticle and quasihole excitation branches are ifTolfTsIl 

i?i,2(p,r) = ^^I^+v(r).p±i?o(p,r), (7) 

where Eq {p, r) = y/[p'^/{2M) - /i(r)]2 + A2(r) is the usual 
spectrum for nonrotating and population balanced mixtures, 
v(r) = ilz X r is the velocity and p is the momentum. In 
Fig.|3] we present three schematic diagrams showing Ei (p, r) 
and i?2(p,r) as a function of p for fixed values of r. At 
each r, the many-body ground state wavefunction fills up all 
of the states with negative energy, and excitations correspond 
to removing a quasiparticle or a quasihole from a filled state 
and adding it to the one that is not filled. In order to show 
that these excitation spectra correspond to three topologically 
distinct superfluid phases that can be observed in atomic sys- 
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tems, next we discuss i?i,2(p, r) in the z = plane such that 

r = {x,y,0). 

First, we consider the population balanced case where /i| = 
In the FPS phase, the excitation spectrum is symmetric 
around the zero energy axis, i.e. Ei{p, r) = —E2{p, r), lead- 
ing to a gapped spectrum as shown in Fig.Oa). However, in 
the rotating case, there is a local asymmetry between the pair- 
ing states {r, p, t; r, -p, 1} and {r, -p, t; r, p, i}. When 
this asymmetry becomes sufficiently large, there exists a mo- 
mentum space region p_ < p < p+ where Ei{—p,r) < 
and E2{p, r) > 0. This occurs for position space region r > 
rx when the condition A{r) = 2Mfl'^r^[^{r) + Mn'^r'^/2] - 
A^(r) > is satisfied, where rx is defined through Airx) = 
and p± = 2M[fi{r) + Am'^r'^] ± 2M^/A{^. Therefore, 
both £'i_2(p, r) have two zeros at p = p±, and the excitation 
spectrum becomes gapless at these momenta. We characterize 
this phase as the PPS, and its excitation spectrum is shown in 
Fig. Kb). 

For a weakly attracting gas, the condition A{r) > can 
only be satisfied near the trap edge when ft <^ uj, but it can 
also be satisfied near the trap center when ft ^ uj. How- 
ever, in the strongly interacting limit where /i is small yet 
positive, this condition can only be satisfied for sufficiently 
small values of A(7') near the trap edge even when ft ^ u. 
Finally, in the molecular limit where fi is negative, this condi- 
tion can not be satisfied anywhere inside the trap, leading to 
a superfluid phase with a gapped excitation spectrum, i.e. the 
Cooper pairs are robust in the molecular limit. As one may ex- 
pect, the asymmetric pairing caused by the rotation does not 
lead to a local population imbalance at any position r since 
Ei{~p,r) < when i?2(p,r) > and vice versa. How- 
ever, this asymmetry prevents the formation of Cooper pairs 
in the phase space region when both i?i 2(p, r) > or both 
i?i.2(p, r) < 0, and thus it is responsible for the creation of 
the PPS and the normal phases. 

'^E,(p,r) / 



\,E,(p,r) 


/ \E,(p,r) 




— 


Eip.r)\ / 


E/p,r)\ 



(a) FPS spectrum 



/ Eip,r) ^ 

(c) PIS spectrum 



FIG. 3: Schematic diagrams showing the excitation spectrum of (a) 
a gapped fully paired superfluid (FPS) phase at r = 0, (b) a gapless 
partially paired superfluid (PPS) phase at r > rx, and (c) a gapless 
population imbalanced superfluid (PIS) phase at r = as a function 
of momentum p. 



We find that the local excitation spectrum changes from 
gapped (FPS) to gapless (PPS) at position r = r^. In homoge- 
nous (infinite) systems such a change is classified as a topo- 
logical quantum phase transition [l^. Since this change oc- 
curs in the momentum space, it does not leave any strong sig- 
nature in the momentum averaged observables such as A(r), 
n(r), etc. However, its direct consequences can be observed 
via a recently developed position and momentum resolved 
spectrocopy 111711 . For instance, the local density distribu- 
tions [jl5ll 



(p, r) = u'ip, r)/[±i?i,2(p, r)] + v\p, r)/[±£;2,i(p, r)], 



are nonanalytic at p — p± when r > ry. Here, u^(p,r) = 
{1 + [e{p) - tiir)]/Eo{p,r)]}/2 and ^{p,r) = {1 - 
[e{p) — /i(r)]/_Eo(p, r)]}/2 are the usual coherence factors 
withe(p) ^p^/{2M). 

For population imbalanced superfluids (PIS) where /i| 7^ 
/ij, the excitation branches shift upwards (downwards), and 
one of them cross zero energy axis when N-j > Ni {N^ < 
Ni), leading to a gapless excitation spectrum. This is ex- 
pected since the excess fermions can only exist in regions with 
both £'i.2(p, r) > or both i?i.2(p, r) < 0. In the absence 
of rotation 1 1511 , for a weakly attracting PIS only one of the 
excitation branches has four zeros as shown in Fig. [3lc) for 
the N'j > Ni case. However, in the strongly attracting limit, 
this branch has only two zeros (not shown). When the sys- 
tem is rotating, the excitation spectra tilt similar to that shown 
in Fig. [3b) (not shown). We remark in passing that a simi- 
lar quantum phase transition with its experimental signatures 
has recently been discussed in the context of trapped p-wave 
superfluids [18]. 

To conclude, we used the BdG formalism to analyze the 
effects of adiabatic rotation on the ground state phases of har- 
monically trapped Fermi gases. We found that the rotation 
breaks Cooper pairs that are located near the trap edge, and 
that this leads to a phase separation between the nonrotating 
superfluid (fully paired) atoms located around the trap cen- 
ter and the rigidly rotating normal (nonpaired) atoms located 
towards the trap edge with a coexistence (partially paired) re- 
gion in between. We also showed that the rotation reveals a 
topological quantum phase transition in the momentum space 
as a function of radial distance. An interesting extention of 
our work is to study emergence of dynamic instabihties for 
fast enough rotation lIsifTgll. 



[1] Ultra-cold Fermi Gases, Procs. of the Int. Sch. of Phys. 'Enrico 
Fermi', Course CLXIV, Varenna 2006, edt. by M. Inguscio, W. 
Ketterle, and C. Salomon (lOS Press, Amsterdam, 2008). 

[2] S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 
80, 1215 (2008). 

[3] N. Nygaard, G. M. Bruun, C. W. Clark, and D. L. Feder, Phys. 



Rev. Lett. 90, 210402 (2003). 
[4] M. Machida and T. Koyama, Phys. Rev. Lett. 94, 140401 
(2005). 

[5] R. Sensarma, M. Randeria, and Tin-Lun Ho, Phys. Rev. Lett. 

96, 090403 (2006). 
[6] M. W. Zwierlein et al. Nature (London) 435, 1047 (2005). 



5 



[7] M. W. Zwieiiein et al. Science 311, 492 (2006). 

[8] I. Bausmerth, A. Recati, and S. Stringaii, Phys. Rev. Lett. 100, 

070401 (2008); and Phys. Rev. A 78, 063603 (2008). 
[9] Hui Zai and Tin-Lun Ho, Phys. Rev. Lett. 97, 180414 (2006). 
[10] Michael Urban and Peter Schuck, Phys. Rev. A 78, 011601(R) 

(2008). 

[11] Y. Ohashi and A. Griffin, Phys. Rev. A 72, 013601 (2005). 
[12] P Castorina et al, Phys. Rev. A 72, 025601 (2005). 
[13] T. Mizushima et al., J. Phys. Soc. Jpn. 76, 104006 (2007). 
[14] L. M. Jensen, J. Kinnunen, and P. Torma, Phys. Rev. A 76, 
033620 (2007). 



[15] M. Iskin and C. A. R. Sa de Melo, Phys. Rev. A 76, 013601 

(2007) . 

[16] G. E. Volovik, cond-mat/0601372 (2006). 

[17] J. T. Stewart, J. P Gaebler, and D. S. Jin, Nature 454, 744 

(2008) . 

[18] M. Iskin and C. J. Williams, Phys. Rev. A 77, 041607(R) 
(2008). 

[19] G. Tonini, F. Werner, and Y. Castin, Eur. Phys. J. D 39, 283 
(2006). 



